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Abstract. We present a new three-dimensional general relativistic hydrodynamics 
code, the Whisky code. This code incorporates the expertise developed over the 
past years in the numerical solution of Einstein equations and of the hydrody- 
namics equations in a curved spacetime, and is the result of a collaboration of 
several European Institutes. We here discuss the ability of the code to carry out 
long-term accurate evolutions of the linear and nonlinear dynamics of isolated 
relativistic stars. 
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1. Introduction 

During the last few years, computa- 
tional general relativistic astrophysics has 
become increasingly important and accu- 
rate. This is partly due to the rapid in- 
crease in computing power through mas- 
sively parallel supercomputers which make 
large-scale, multi-dimensional numerical 



simulations possible ( 


Font et al. 


2000b 


Shibata fc Urvu 


20001: iFont et al. 


2002 


Shibata & Urvu 


2002; LDuez et al 2002) 



In addition to being a unique tool for in- 
vestigating General Relativity in regimes 
of strong and rapidly varying gravita- 
tional fields, such simulations are also 
needed to fully understand the incoming 
wealth of observations from high-energy 
astronomy and (near-future) gravitational 
wave astronomy. In an attempt to re- 
spond, at least in part, to this need 



we have recently developed Whisky, a 
three-dimensional (3D) general relativis- 
tic hydrodynamics code. The Whisky code 
is the result of an ongoing collabora- 
tion among several European Institutes, 
i.e. the Albert Einstein Institute (Golm, 
Germany), SISSA (Trieste, Italy) and the 
Universities of Thessaloniki (Greece) and 
Valencia (Spain), joined in a European 
Network investigating sources of gravita- 
tional waves (see www . eu-network . org for 
further information). 

In practice, the Whisky code solves the 
general relativistic hydrodynamics equa- 
tions on a 3D numerical grid with 
Cartesian coordinates. The code has been 
constructed within the framework of 
the Cactus Computational Toolkit (see 
www.cactuscode.org for details), devel- 
oped at the Albert Einstein Institute 
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(Golm), which provides high-level facilities 
such as parallelization, input/output, etc., 
but also the solution of the Einstein equa- 
tions with matter terms being provided by 
the Whisky code. 

In this paper we briefly discuss the main 
features of Whisky and present some tests 
of its validation. Indeed, validation repre- 
sents an important aspect of the develop- 
ment of any modern 3D finite-difference 
code. The reasons for this are rather sim- 
ple and are related to: (i) the lack of pre- 
cise knowledge of the space of solutions of 
the coupled system of the Einstein and gen- 
eral relativistic hydrodynamics equations; 
(ii) the likely chance that coding errors are 
made in the implementation of the thou- 
sands of terms involved in the solution of 
such a complicated set of coupled partial 
differential equations; (Hi) the complexity 
of the computational infrastructure needed 
for the use of the code in a massively par- 
allel environment which increases the risk 
of computational errors. 

The tests presented here will show 
both the accuracy and the convergence 
of our formulation for the general rel- 
ativistic hydrodynamics equations, which 
are coupled to a conformal transverse- 
tracele ss formulation of the Ei nstein equa- 
tions (jNakamura et al. I 119871 ). They will 



also show the ability of the code to follow 
stably the linear and nonlinear dynamics 
of isolated relativistic stars. More specif- 
ically, we will first present results of the 
linear pulsations of spherical and rapidly 
rotating stars. The computed frequencies 
of radial and quasi-radial oscillations will 
be compared with the corresponding fre- 
quencies obtained with lower-dimensional 
numerical codes, or with alternative tech- 
niques such as the Cowling approximation 
(in which the spacetime is held fixed and 
only the general relativistic hydrodynam- 
ics equations are evolved), or with relativis- 
tic perturbative methods. The comparison 
shows an excellent agreement confirming 
the ability of the code to extract physically 
relevant information even from tiny per- 
turbations. The successful determination 



of the eigenfrequencies for rapidly rotating 
stars computed with our code is particu- 
larly noteworthy since the equivalent prob- 
lem has not yet been tackled with pertur- 
bative techniques. 

We will also investigate the nonlinear 
dynamics of stellar models that are un- 
stable to the fundamental radial mode of 
pulsation. We show that upon perturba- 
tion, the unstable models will either col- 
lapse to a black hole, or migrate to a con- 
figuration in the stable branch of equilib- 
rium configurations. In the case of a grav- 
itational collapse, we will follow the evo- 
lution all the way down to the formation 
of a black hole, tracking the generation of 
its apparent horizon. In the case of mi- 
gration to the stable branch, on the other 
hand, we will be able to accurately fol- 
low the nonlinear oscillations that accom- 
pany this process and that can give rise 
to strong shocks. The ability to simulate 
large amplitude oscillations is important as 
we expect a neutron st ar formed in a su- 



pernova core-collapse (|Zwerger fc Miiller . 
1997 : Dimmelmeier et al. 200lh or in the 



accretion-induced collapse of a white dwarf 
to oscillate violently in its early stages of 
life. 

We use the signature (—,+,+,+) and 
units in which c = G = M@ = 1. Greek in- 
dices are taken to run from to 3 and Latin 
indices from 1 to 3. 



2. Basic Equations 

The Whisky code has been constructed 
exploiting the expertise developed in 
the building of a similar but distinct 
code whi ch has been desc r ibed exten- 
sively in (front et al. I l2000bl |2002h . The 
development of the Whisky code has 
benefited from the release of a public 
version of the general relativistic hy- 
drodynamics code described in Font 
et al. (2000b), (2002), and developed 
mostly by the group at the Washington 
University as part of the NASA Neutron 
Star Grand Challenge Project (see 
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wugrav.wustl.edu/Codes/GR3D for de- 
tails). 

The Whisky code, however, also 
incorporates important recent devel- 
opments regarding, in particular, new 
numerical methods for the solution of 
the hydrodynamics equations. These 
include: (i) the Piecewise Parabolic 
Method (PPM) dColella fc Woodward 



1984[ ) and the Essen tially Non- Osc il latory 



(ENO) methods (|Harten et al. I 119871) 



for the cell reconstruction procedure; 
(ii) the Harten-Lax-van Leer-Einfcldt 



(HLLE) (jrlarten et al. I Il983l) approxi 



mate Rie mann solver , the M arquina flux 
formula (jAlov et al. I Il999bl ); (Hi) the 
analytic expression fo r the left eigenvectors 



( Ibafiez et al. _ 2000l) and c ompact flux 



formulae (jAlov et al. I Il999al ) for the Roe 
Riemann solver and the one using the 
Marquina flux formula; (iv) the possibility 
to couple the general relativistic hydro- 
dynamics equations with a conformally 
decomposed 3-metric. The incorporation 
of these new numerical techniques in the 
code has lead to a much improved ability 
to simulate relativistic stars, as will be 
shown in the section devoted to the tests. 
The interested reader m ay also refer to 
Font et al. I (|2000bl . l2002h for more details 
about the general formulation and the 
results obtained with the similar, distinct 
code. 

While the Cactus code provides at each 
time step a solution of the Einste in equa- 
tions (see Alcubierre et al. ( 2000l) for the 
validation of the Cactus code for the space- 
time evolution) 



Giiu — 87rT, lt/ 



(1) 



where G^ v is the Einstein tensor and T^ v is 
the stress-energy tensor, the Whisky code 
provides the time evolution of the hy- 
drodynamics equations, expressed through 
the conservation equations for the stress- 
energy tensor T^ v and for the matter cur- 
rent density J M 



An important feature of the Whisky 
code is the implementation of the so 
called Valencia formu lation of the hy- 
drodynamics equatio ns dMart i et al. Ill99ll 



Banvuls et ail Il997i: llbanez et al. I l200dh 



in which the set of equations ([2]) is writ- 
ten in hyperbolic, flux-conservative form 



&q + a,fW(q)=B(q), 



(3) 



where the right hand side (the source 
terms) depends only on the metric and 
on the stress-energy tensor, vanishing in a 
flat spacetime, where the strict hyperbol- 
icity is recovered. In order to write system 
([2]) in the form of system ((3]), the 'prim- 
itive hydrodynamical variables (i.e. the 
rest mass-density p, the pressure p, the 
fluid 3-velocities v l , the internal energy 
density e and the Lorentz factor W) are 
mapped to the so called conserved variables 
q = (D,S 1 ,t) via the relations 



D = ^Wp , 
S l = ^ P hW 2 v l , 
r = V7 (phW 2 -p)-D 



(4) 



where 7 is the determinant of the spatial 3- 
metric 7^ and h = 1 + e+p/ p is the specific 
enthalpy. The Lorentz factor is defined in 
terms of the velocities and of the 3-metric 
as W = (1 - 7i 3 vV)~ 1/2 . Note that only 
five of the primitive variables are indepen- 
dent. 

Finally, an equation of state is used to 
relate the pressure to the rest mass density 
and to the energy density. The code can use 
any equation of state, but presently only a 
polytropic equation of state of the type 

P = Kp T , (5) 
an "ideal fluid" equation of state 



P=(r-l)pe 



(6) 



a nd a hybrid equa tion of state as described 



in lZwerger I (jl995l ) have been implemented. 
For more details on this formulation, see 







V M J" = 0. 



(2) 



also the review bv lFontl (l2000h . 

An important feature of the Valencia 
formulation is that it is possible to extend 
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to relativistic hydrodynamics the power- 
ful numerical methods developed in classi- 
cal hydrodynamics; in particular, our code 
takes advantage of the High Resolution 
Shock Capturing (HRS C) properties o f 



Godunov type methods ( Godunov 1959h 



For a f ull intr o ductio n to HRSC m eth- 
ods see lLanev I (|l998l ). iToro I (|l999h and 



LeVeaue I (1998) 



3. Numerical Techniques 

The time update of all the equa- 
tions, general relativistic hydrodynamics 
and Einstein, are pe rformed with the 
Meth od of Lines (MoL) (jLanev Hl998l:lToro 



1999). The method of lines is a procedure to 
transform a set of partial differential equa- 
tions such as Q into a set of ordinary dif- 
ferential equations. This is done by inte- 
grating equations ([3]) over space in every 
computational cell defined by its position 
(xi,Vj,z k ) 



:(q) = L(q) 



dt 

Vj + l/2 f z k + l/2 



sd 4 x 



f (1) (q(x i _i /2 ,y,z))dydz 



Vj-1/2 J z k~l/2 
Vj + l/2 f z k + l/2 



rVj+i/2 rZk+1/2 

f { >(q(x l+1/2 ,y,z))dydz 

J Vj-1/2^ z k-l/2 



+ 



(7) 



where q is the spatially integrated vector 
of conserved variables, i.e. 



q da; dy dz 



(8) 



and f W is the i-th component of the flux 
five- vector f . 

Several time- integrators are available in 
our implementation of MoL and the or- 
der of accuracy of the solution of the ordi- 
nary differential equation (|7J| is the same as 
the truncation order of the integrator em- 
ployed, provided that the discrete operator 
L is of the same order in space and at least 
first-order accurate in time. 

In our implementation of MoL, the right 
hand side operator L(q) is simplified by 



approximating the integrals with the mid- 
point rule to get 



L(q) — Sij,k + ff_i/ 2 j,fe *i+i/2,j,fc + • ■ • ■ 

(9) 

Given this simplification, the calculation of 
the right hand side of (JT]) splits into the 
following parts: 



(i) 



1. 



Calculation of the source terms 
s(q(xj, t/j, Zk)) at all the grid points. 
Reconstruction of the data q to both 
sides of a cell boundary. In this way, two 
values q^ and q R of qi+i/2,j,fe are deter- 
mined at cell boundary; q^ is obtained 
from cell i (left cell) and q R from cell 
i + 1 (right cell) . The code implements 
several reconstruction methods. In par- 
ticular, as Total Variation Diminishing 
(TVD) methods we have implemented 
"minm od" , van Leer m onotonized cen- 
tered (y^mLeer. 19791 ) and Superbee 



(jToro I 119991 ). Additional reconstruc- 
tion methods are : arbitrary order 



ENO methods (jHarten et al. I 119871 ) 



and the Piecewise Parabolic Method 



(|Colella fc Woodward 1119841 ). which is 
a third order accurate in space. As men- 
tioned below, PPM has emerged as our 
actual best choice for all the test evolu- 
tions we present here. 
Solution at cell b oundary of the 



Riemann problem (ILeVeaue 1998; 



Torol Il999t lLanev I Il998l ) having the 



values q^ R as initial data. 

4. Calculation in each coordinate di- 
rection (x, y, z) of the inter-cell 
flux f (a:) (q l+ i/2j,fe), f (y) (qij+i/2,fe), 
f ' < " z H c li.j,k+i/2) 7 that is the flux across 
the boundary between a cell (e.g. the 
i-th) and its closest neighbour (e.g. the 
(* + l)-th). 

5. Recovery of the primitive variables and 
computation of the stress-energy tensor 
for use in the Einstein equations. 

As a result of steps 1. - 4., the core of 
the Whisky code is effectively represented 
by two routines. One that reconstructs the 
function q at the boundaries of a computa- 
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tional cell and another one that calculates 
the inter-cell flux f at this cell boundary. 

As for the reconstruction methods, a 
similar variety is present for the approx- 
imate Riemann solvers implemented in 
the Whisky code. More s pecifically, are 
availa ble the fast HLLE (H arten et al 



19831) and the widespread Roe (|Roe 1981 ) 
solvers as well as the accurate M arquina 
flux formula (|Alov et al. I Il999bl) . which 



is used to solve the Riemann problem in 
a way that differs from the Roe solver 
only at sonic points, where the Roe solver 
has problems. The Roe-based approximate 
Riemann solvers need the computation 
of the eigenvalues and eigenvectors (from 
both the right and left cell) of the linearized 
Jacobian matrices A L and A H given by 
i L = A L q L and = A R q R . We have im- 
plemented the a nalytic expression f or the 
left eigenvectors ( Ibahez et al. 2000), thus 
avoiding the computationally expensive in- 
version of the three 5x5 matrices of the 
right eigenvectors, associated to each spa- 
tial direction. This implementation brings 
a 40% reduction of the computational time 
spent in the solution of the hydrodynamics 
equations. However, in evolutions involving 
also the time integration of the Einstein 
equations, this is reduced to a 5% decrease 
in computational cost. This is due to the 
fact that the largest part of the time is 
spent in the update of the spacetime field 
variables. 



4. Numerical Tests 

As mentioned above several tests have 
been performed to assess the stability and 
accuracy of our code. The results obtained 
so far are very encouraging and already 
during these preliminary steps some new 
physical results have been achieved. 

First of all, we consider a standard 
shock tube test, setting as initial data a 
global Riemann problem, i.e. one in which 
the initial discontinuity is orthogonal to the 
main diagonal of the cubic grid. More pre- 
cisely the initial data consist of a left and 




Fig. 1. Solution of a Riemann problem set 
on the main diagonal of the cubic grid. The 
figure shows the comparison of the hydro- 
dynamical variables evolved by Whisky, in- 
dicated with symbols, with the exact so- 
lution. The numerical simulation was ob- 
tained with the van Leer reconstruction 
method and the Roe solver, on a 140 3 grid. 

right state with values 

PR = 1; PR = 1-666 x 10~ 6 ; v R = 



p L =10; pl = 1.333; 



v L =0 



In Fig. Q] we show the solution at a 
given time together with the exact solu- 
tion. The excellent agreement of the two 
sets of curves is particularly remarkable if 
one bears in mind that the initial shock 
is placed on the main diagonal of the cu- 
bic grid, so that the evolution is fully 3- 
dimensional. 

Next, we consider the evolution of sta- 
ble relativistic polytropic spherical (TOV) 
stars. As this is a static solution, no evolu- 
tion is expected. Yet as shown in Fig. [3J 
both a small periodic oscillation and a 
small secular increase of the central density 
of the star are detected during the numeri- 
cal evolution of the equations. Both effects 
have, however, a single explanation. Since 
the initial data contains also a small trun- 
cation error, this is responsible for trigger- 
ing radial oscillations which appear as pe- 
riodic variations in the central density. As 
the resolution is increased, the truncation 
error is reduced and so is the amplitude 
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of the oscillation. The secular growth, on 
the other hand, is a purely numerical prob- 
lem, probably related to the violation of 
the constraint equations. As for the oscilla- 
tions, also the secular growth converges to 
zero with increasing resolution. The conver- 
gence properties of the code are also clearly 
shown in the growth of the Hamiltonian 
constraint violation (Fig. [3]) , where we can 
see that almost second-order convergence is 
achieved. Note that the convergence rate is 
not exactly second-order, because the re- 
constructi on schemes are only first-order 
accurate ( Alcubierre et al. 1 120001) at local 
extrema (i.e. the center and the surface of 
the star) thus increasing the overall trun- 
cation error. 

As anticipated, the PPM reconstruction 
scheme has shown to be more accurate than 
the TVD ones. This is clear in Fig. [H in 
which the results obtained with the PPM 
method are compared to the best of the 
TVD methods (i.e. the van Leer one) for 
a stable TOV run, similar to the previous 
ones, using 64 3 grid points. Note that the 
PPM reconstruction is more effective in re- 
ducing both the initial truncation error (as 
shown by the smaller amplitudes in the os- 
cillations) and the secular error (as shown 
by the small growth rate). 

In order to further investigate the accu- 
racy of our implementation of the hydrody- 
namics equations, we have suppressed the 
spacetime evolution and solved just the hy- 
drodynamics equations in the fixed space- 
time of the initial TOV solution. This ap- 
proximation is referred to as the "Cowling 
approximation" and is widely used in per- 
turbative studies of oscillating stars. In 
this case, in addition to the confirmation 
of the convergence rate already checked 
in fully evolved runs, we have also com- 
pared the frequency spectrum of the nu- 
merically induced oscillation with the re- 
sults obtained by a n independent 2D code 
(|Font et al. Il2000al) and with perturbative 
analyses. In Fig. [5] we show a compari- 
son between the two codes reporting the 
power spectrum of the central density os- 
cillations computed with the Whisky code 



32 




64 




96 




128 

/ ' 













t (ms) 



Fig. 3. The L2-norm of the Hamiltonian 
constraint for the same evolutions as in 
Fig.H 



and the corresponding frequencies as ob- 
tained with perturbative techniques and 
with the 2D code. Clearly the agreement 
is very good with an error below 1% in the 
fundamental frequency. The fact that the 
frequencies computed with the code coin- 
cide with the physical eigenfrequencies cal- 
culated through perturbative analysis al- 
lows us to study with our code the physical 
properties of linear normal-modes of oscil- 
lation even if such oscillations are generated 
numerically. 

The last test performed in the linear 
regime consisted of the evolution of rapidly 
rotating stars, with angular velocity up to 
95% of the allowed mass-shedding limit for 
uniformly rotating stars. The initial data 
routin es have been adapted from th e RNS 
code (jStergioulas fc Friedman 1 Il995h . As 
in the previous tests the Hamiltonian con- 
straint shows a convergence rate of nearly 
second-order everywhere, except at the sur- 
face and the center of the star. In analogy 
with the nonrotating case, the truncation 
error triggers quasi-radial oscillations in the 
star. Such pulsations converge to zero with 
increasing resolution. Determining the fre- 
quency spectrum of fully relativistic and 
rapidly rotating stars is an important sci- 
entific achievement, allowing the investiga- 
tion of a parameter space which is astro- 
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Fig. 2. Central mass-density, normalized to the initial value, in a stable TOV star 
(M = 1.4M Q and polytropic index T = 2) evolution at different resolutions. PPM and 
Marquina were used for all runs. 



physically relevant but too difficult to treat 
with current perturbative techniques. 

Note that a number of small improve- 
ments on the boundary and gauge condi- 
tions have allowed us to extend consider- 
ably the timescale of our evolutions of sta- 
ble rapidly rotating stars, which can now 
be evolved for about 10 ms, that is for sey - 



eral rotational periods (St ergioulas II2002T) 



We now consider tests of the nonlinear 
dynamics of isolated spherical relativistic 
stars. To this purpose we have constructed 
TOV solutions that are placed on the un- 
stable branch of the equilibrium configu- 
rations (see inset of Fig. [5]). The trunca- 
tion error in the initial data for a TOV is 
sufficient to move the model to a different 



configuration and in Whisky this leads to 
a rapid migration toward a stable configu- 
ration of equal rest-mass but smaller cen- 
tral density. Such a violent expansion pro- 
duces large amplitude radial oscillations in 
the star that are either at constant ampli- 
tude if the polytropic equation of state §5§ 
is used, or are damped through shock heat- 
ing if the ideal fluid equation of state ^ is 
used. This is summarized in Fig. [51 which 
shows the time series of the normalized cen- 
tral density for a TOV. We also show that 
the asymptotic central density tends to a 
value corresponding to a rest-mass slightly 
smaller than the initial one (straight line). 
This is the energy loss due to the internal 
dissipation. 
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12 3 

t (ms) 

Fig. 6. Normalized central mass-density evolution of an M = 1.4 Mq, T — 2 unstable 
TOV star performed with 96 3 grid points. 



An alternative solution for the unstable 
model is that of a gravitational collapse. 
In this case, in fact, the initial star does 
not expand but rather moves to increas- 
ingly larger density configurations, finally 
forming a black hole. In order to study 
the gravitational collapse of the unstable 
configuration, the introduction of a density 
perturbation in the initial model is nec- 
essary. A very small one of the order of 
1% with dependence cos(7rr/2r s ), where r 
is the coordinate distance from the center 
and r s its value on the surface, is suffi- 
cient to overcome the effects of the trunca- 
tion error and induce the star to collapse. 
Note that after adding the perturbation 
to the initial configuration, the constraint 
equations are solved to provide initial data 



which are a solution to the Einstein equa- 
tions. As a summary of the results ob- 
tained, we show in Fig. [7] the growth of 
the horizon mass, tracked with an apparent 
horizo n finder based on the fast-flow algo- 
rithm (jGundlach Ill998h . At t = 0.246 ms a 



black hole forms and an apparent horizon 
appears. As the remaining stellar material 
continues to accrete onto the newly formed 
black hole, its horizon mass increases, fi- 
nally leveling off until t — 0.306 ms. The 
subsequent growth of the horizon mass is 
then only the result of the increasing er- 
ror due to grid stretching (the radial met- 
ric function develops a sharp peak which 
cannot be resolved accurately enough). 
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t(ms) 

Fig. 4. Central mass-density, normalized 
to the initial value, in a stable TOV star 
(M = 1.4 M© and polytropic index T = 2) 
evolution with 64 3 grid points. Comparison 
between PPM and van Leer reconstruction 
methods. 



0.30 
0.25 
0.20 
0.15 
0.10 
0.05 




Fig. 5. Fourier transform of the central 
mass-density evolution of an M — 1.4 Mq, 
T = 2 stable TOV star performed with 128 3 
grid points. The units of the vertical axis 
are arbitrary. 



5. Conclusions 

We have illustrated the main features 
and the present status of our new 3D 
general relativistic hydrodynamics code. 
Through a wide set of numerical tests, we 
have shown both the accuracy and the con- 
vergence of our implementation of the for- 
mulation of the general relativistic hydro- 
dynamics equations, coupled to a confor- 
mal transverse-traceless formulation of the 



1.3 ' — 
0.23 



0.27 
t(ms) 



Fig. 7. Horizon mass evolution. The initial 
rest mass of the TOV star is M = 1.44 M© 
and again T = 2. 



Einstein equations. We have also shown 
that our code can accurately and stably 
evolve the linear and nonlinear dynamics 
of isolated relativistic stars, both for pul- 
sations of spherical and rapidly rotating 
stars. The computed frequencies of radial 
oscillations are compared with the corre- 
sponding frequencies obtained with other 
numerical and perturbative techniques and 
the good agreement among these values is 
one of the greatest present achievements of 
our code. 

We have also investigated the nonlinear 
dynamics of stellar models that are unsta- 
ble to the fundamental radial mode of pul- 
sation, decaying either to a black hole (a 
collapse which we can follow up to and past 
the formation of the event horizon) or in a 
migration to a configuration on the stable 
branch of equilibrium configurations. 

These encouraging results are impor- 
tant premises for the application of the 
Whisky code to more physical scenarios. 
Furthermore we are presently working on 
incorporating in Whisky the use of a fixed 
mesh refinement. This is an important im- 
provement that will increase the numerical 
resolution where needed and move the po- 
sition of the outer boundaries further out 
in the wave zone, where information on the 
gravitational wave content of the spacetime 
can be reliably extracted. 
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